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Abstract 

For  more  than  a  decade  the  US  government  has  been  developing  laser-based  sensors  for 
detecting,  locating,  and  classifying  aerosols  in  the  atmosphere  at  safe  standoff  ranges. 
The  motivation  for  this  work  is  the  need  to  discriminate  aerosols  of  biological  origin  from 
interferent  materials  such  as  smoke  and  dust  using  the  backscatter  from  multiple 
wavelengths  in  the  long-range  IR  (LWIR)  spectral  region.  Through  previous  work, 
algorithms  have  been  developed  for  estimating  the  aerosol  spectral  dependence  and 
concentration  range-dependence  from  these  data.  The  range-dependence  is  required  for 
locating  and  tracking  the  aerosol  plumes,  and  the  backscatter  spectral  dependence  is  used 
for  discrimination  by  a  support  vector  machine  classifier.  Substantial  progress  has  been 
made  in  these  algorithms  for  the  case  of  a  single  aerosol  present  in  the  lidar  line-of-sight 
(LOS). 

Often,  however,  mixtures  of  aerosols  are  present  along  the  same  LOS  overlapped  in 
range  and  time.  Analysis  of  these  mixtures  of  aerosols  presents  a  difficult  inverse 
problem  that  cannot  be  successfully  treated  by  the  methods  used  for  single  aerosols. 
Fortunately,  recent  advances  have  been  made  in  the  analysis  of  inverse  problems  using 
shrinkage-based  Li -regularization  techniques.  Of  the  several  Li -regularization  methods 
currently  known,  the  split  Bregman  algorithm  is  straightforward  to  implement,  converges 
rapidly,  and  is  applicable  to  a  broad  range  on  inverse  problems  including  our  aerosol 
unmixing.  In  this  paper  we  show  how  the  split  Bregman  algorithm  can  successfully 
resolve  LWIR  lidar  data  containing  mixtures  of  bioaerosol  simulants  and  interferents  into 
their  separate  components.  The  individual  components  then  can  be  classified  as  bio-  or 
non-bio  aerosol  by  our  SVM  classifier.  We  illustrate  the  approach  through  data  collected 
in  field  tests  over  the  past  several  years  using  the  US  Anny  FAL  sensor  in  testing  at 
Dugway  Proving  Ground,  UT. 
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1.  Introduction 


This  work  describes  a  generalization  of  previous  research1 2  on  developing  efficient 
algorithms  for  estimating  the  parameters  of  optically  thin  aerosols  in  the  atmosphere 
using  data  from  rapidly  tuned  multiple-wavelength  long  wave  infrared  (LWIR)  lidar. 
The  motivation  for  this  work  remains  the  same:  the  need  for  detecting,  locating,  tracking, 
and  discriminating  atmospheric  aerosols  at  safe  standoff  ranges  using  time-series  data 
collected  at  a  discrete  set  of  CO?  laser  wavelengths.  Our  earlier  work  considered  a  single 
aerosol  material  producing  a  backscatter  enhancement  over  that  from  the  natural 
atmosphere.  The  goals  were  to  detect  and  track  the  aerosol  by  means  of  estimates  of  the 
range-dependence  of  the  concentration,  and  to  discriminate  potentially  harmful 
aerosols — particularly  those  of  biological  origin — from  interferents  such  as  smoke  and 
dust  by  means  of  the  spectral  dependence  of  the  backscatter.  State-of-the-art  support 
vector  machine  (SVM)  classifiers  have  been  developed  for  this  discrimination. 

Often,  however,  mixtures  of  aerosols  are  present  along  the  same  lidar  line-of-sight 
overlapped  in  range  and  time.  For  example,  it  is  quite  possible  to  have  a  biological  or 
chemical  agent  release  from  a  munition  accompanied  by  dust  and  other  byproducts  of  the 
explosive  release.  Those  extra  materials  can  distort  the  spectral  dependence  of  the  threat 
material  and  degrade  the  ability  of  the  sensor  to  correctly  identify  the  threat.  The  analysis 
of  aerosol  mixtures  presents  a  difficult  inverse  problem  that  cannot  be  successfully 
treated  by  the  methods  used  previously  for  single  materials. 

Fortunately,  recent  advances  have  been  made  in  the  analysis  of  inverse  problems  using 
shrinkage-based  L 1 -regularization  techniques.  We  focus  here  on  the  split  Bregman 
algorithm"  because  it  is  straightforward  to  implement,  converges  rapidly,  and  is 
applicable  to  a  broad  range  of  inverse  problems  including  our  aerosol  unmixing.  In  this 
paper  we  show  how  the  split  Bregman  algorithm  can  successfully  resolve  LWIR  lidar 
data  containing  mixtures  of  bioaerosol  simulants  and  interferents  into  their  separate 
components.  The  individual  components  then  can  be  classified  as  bio-  or  non-bio  aerosol 
by  our  SVM  classifier. 

The  use  of  Bregman  iteration  to  solve  difficult  Li -estimation  problems  efficiently  was 
described  by  Yin  et  a  1. 3  In  that  paper  they  showed  the  equivalence  of  Bregman  iteration 
to  the  augmented  Lagrangian  method  for  the  minimization  of  convex  functions  with 
linear  constraints.  In  a  subsequent  paper,  Goldstein  and  Osher  [2]  showed  that  Bregman 
iteration  could  be  efficiently  implemented  by  splitting  the  original  combination  of  Li  and 
L2  optimizations  over  a  single  variable  into  two  separate  minimizations;  one  over  a 
differentiable  quadratic  component,  and  a  second  over  a  non-differentiable  Li-norm 


1  R.  E.  Warren,  R.  G.  Vanderbeek,  A.  Ben-David,  and  J.  L.  Ahl,  Simultaneous  estimation  of  aerosol  cloud 
concentration  and  spectral  backscatter  from  multiple-wavelength  lidar  data,  Appl.  Optics,  47,  pp  4309- 
4320,  2008. 

2  T.  Goldstein  and  S.  Osher,  The  split  Bregman  method  for  Ll-regidarized problems.  SIAM  J.  Imaging  Sci., 
2(2):323-343,  2009.  ISSN  1936-4954. 

3  W.  Yin,  S.  Osher,  D.  Goldfarb,  and  J.  Darbon,  Bregman  iterative  algorithms  for  I /-minimization  with 
application  to  compressed  sensing,  SIAM  J.  Imaging  Sci.,  vol.  1,  pp  143-168,  2008. 
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component  with  an  objective  function  that  could  be  solved  in  closed  form  by  shrinkage. 
This  is  the  origin  of  the  term  split  Bregman.  Subsequently,4’5  it  was  recognized  that  their 
splitting  algorithm  could  be  interpreted  from  a  Lagrangian  and  penalty  standpoint 
(augmented  Lagrangian)  analogous  to  that  described  in  [3].  It  is  this  augmented 
Lagrangian  interpretation  of  split  Bregman  that  we  follow  here.  For  a  comprehensive 
analysis  of  the  interrelations  between  Bregman  iteration,  augmented  Lagrangians,  split 
Bregman,  and  other  techniques,  see  the  review  by  Esser.6 

In  our  treatment  of  single  aerosol  estimation,  we  gave  a  rather  detailed  discussion  of  the 
lidar  measurement  and  preprocessing  needed  to  remove  the  effects  of  the  natural 
atmosphere  and  to  deconvolve  the  long  CO?  transmitter  pulse  waveforms.  Those 
discussions  pertain  here  without  any  modification.  We  therefore  concentrate  on  the 
problem  of  generalizing  the  estimation  to  multiple  aerosols  by  use  of  the  split  Bregman 
method.  Accordingly,  section  2  provides  a  brief  summary  of  the  advantages  of  Li- 
regularization  over  the  more  traditional  L2  method,  and  more  particularly,  how  this 
newer,  more  robust  method  can  be  implemented  efficiently  through  split  Bregman. 
Section  3  describes  how  the  multi-aerosol  estimation  is  fonnulated  as  a  dual  application 
of  the  basic  split  Bregman  algorithm.  In  section  4  we  give  some  examples  of  processing 
aerosol  mixture  data  from  both  simulated  and  actual  release  data  collected  by  the  US 
Army  FAL  sensor  during  testing  at  Dugway  Proving  Ground,  UT.  Section  5  summarizes 
and  concludes. 

2.  Li-Regularization  of  Inverse  Problems  by  the  Split  Bregman  Method 

It  is  well  known  that  inverse  problems  such  as  the  retrieval  of  signals  observed  after  some 
smoothing  operation  such  as  blurring  are  ill-posed  in  the  sense  that  small  changes  in  the 
data  due  to  noise,  for  example,  can  have  large  effects  on  the  solution.  Otherwise  stated, 
many  different  solutions  are  compatible  with  the  same  data.  Since  most  of  these 
“solutions”  are  wildly  oscillatory,  and  therefore  not  physically  meaningful,  some  sort  of 
regularization  is  required  to  produce  useful  results.  Historically,  L2-regularization  has 
been  used  in  the  context  of  least-squares  estimation  since  it  is  easy  to  implement  through 
ridge  regression.  Although  stabilizing  the  solution,  ridge  regression  tends  to  be 
unselective  in  the  features  it  suppresses  and  can  lead  to  severely  biased  results. 

As  noted  above,  more  recent  alternatives  to  L2 -regularization  have  replaced  the  quadratic 
regularization  term  in  the  objective  function  with  an  Li-norm  regularizing  term.  This  has 
the  effect  of  selectively  removing  features  in  the  estimate  that  produce  the  large 
fluctuations  associated  with  direct  inversion  without  regularization.  This  selectivity  gives 
L  i  estimates  a  sparse  representation  in  an  appropriate  basis  set.  The  downside  of  L  i  - 


4  X.  Tai  and  C.  Wu,  Augmented  Lagrangian  method,  dual  methods  and  split  Bregman  iteration  for  ROF 
model.  Scale  Space  and  Variational  Methods  in  Computer  Vision:  Lecture  Notes  in  Computer  Science,  vol. 
5567/2009,  pp  502-513,  2009. 

5  Y.  Gui-Bo  and  X.  Xiaohui,  Split  Bregman  for  large  scale  fused  Lasso,  J.  Computational  Statistics  and 
Data  Analysis,  vol.  55,  no.  4,  April,  2011. 

6  E.  Esser,  Applications  of  Lagrangian-based  alternating  direction  methods  and  connections  to  split 
Bregman,  UCLA  CAM  report  09-31,  2009. 
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regularization  has  been  the  loss  of  differentiability  of  the  objective  function  making 
implementation  more  difficult  than  L2. 

As  noted  above,  the  new  class  of  shrinkage-based  inversion  methods  such  as  split 
Bregman  largely  overcomes  the  computational  limitations  of  Li -regularization.  To 
illustrate  the  advantages  of  Li  over  L2  we  compare  them  on  a  simple  one-dimensional 
simulated  example.  Our  objective  is  to  minimize  the  functionals  Jx(u)  and  J2(u) : 

«,  2  =  argmin  Jl2{u) , 

’  u  ’ 

1  ..  ||2  ,,  I. 

Jx(u)  =  —  jKu  ~  f\\2  +  A  ||w||j  5 

where  J\^  are  the  sum  of  a  quadratic  term  quantifying  the  fidelity  of  the  linear  model  Ku 
to  the  data  /  and  the  second  terms  are  regularizers  using  the  Li  and  L2  norms, 
respectively: 


1 =11 


mL  = 


II 


-11/2 


We  assume  the  model  operator  K  and  the  data  to  be  known.  The  scalar  coefficients  X  and 
//  are  set  heuristically  to  balance  the  conflicting  needs  for  fitting  the  model  and  solution 
smoothness. 


For  concreteness  we  modeled  u  as  the  sum  of  two  narrow  Gaussian  functions  given  by 


w(x)  =  exp  -(jc  +  O.l)”  / 0.001  +exp  -(x-0.l)2  / 0.00 1 


and  shown  in  Figure  1  sampled  uniformly  over  1000  points  between  -2  and  2.  The  data  / 
were  created  by  blurring  u  with  a  Toeplitz  matrix 


with  r  =  0.99,  and  adding  zero-mean  independently  distributed  (pseudo-)random  normal 
noise  with  standard  error  a  =  0.004: 

■fn  ^  '  1  ^ ^  9 n  1 

m 

q n  D  ^(0,a2). 
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The  simulated  data  are  shown  in  Figure  1  together  with  the  L2  regression  estimate  at  the 
value  //  =  0.015,  and  the  Li  estimate  computed  by  the  split  Bregman  method  defined 
below  also  using  A  =0.015.  We  note  the  L2  estimate  is  unable  to  resolve  the  peaks, 
whereas  the  Li  estimate  gives  a  clean  separation  of  the  two  components. 


Fig.  1  Data  and  estimates  from  the  narrow  Gaussian  example. 


The  objective  function  J\  is  non-smooth  because  of  the  Li-nonn  regularization  term. 
This  makes  the  optimization  over  u  difficult  by  the  usual  calculus-based  methods.  The 
basic  idea  of  split  Bregman  is  to  introduce  additional  variables  into  the  J\  objective 
function  that  will  allow  a  convenient  separation  of  the  minimizations  over  the  quadratic 
(smooth)  data  tenn 


H(u)^\Ku-f\l 

and  the  non-smooth  regularization  term  ||w||  .  This  is  done  in  two  steps.  First,  the 
unconstrained  problem  arg  min  II(ii)  +  M  is  replaced  by  the  equivalent  constrained 
problem 

arg  min  H (w)  +  ||d||  such  that  u  =  d  . 

1 1  ^  m 


Second,  the  constraint  u  =  d  is  enforced  by  adding  a  Lagrange  multiplier  term  and  an 
additional  quadratic  penalty  term  to  get  the  Lagrangian 


L[u,d,b^  =  +  +  A  (b,u—d}  + 


(1) 
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where  b  is  the  Lagrange  multiplier  vector,  (a,b)  denotes  the  inner  product  of  vectors  a 

and  b,  and  A  is  a  positive  regularization  parameter.  In  this  formulation,  L  assumes  the 
form  of  an  augmented  Lagrangian  [4,  5]  to  be  minimized  for  u  and  d,  and  maximized  for 
the  dual  variable  b.  In  the  appendix  the  derivation  of  this  saddle  point  optimization  is 
sketched.  The  resulting  algorithm  is  given  as  Algorithm  1.  The  function  SA  (x) 

appearing  in  Algorithm  1  is  the  shrinkage  operator  that  plays  a  key  role  in  the  success  of 
the  method  by  promoting  sparseness  in  the  solution  through  shrinking  input  values 
toward  zero.  S  is  defined  as 


x  —  X,x  >  A 

SA(x)  =  < 

0,-X  <  x  <  A 

>  . 

x  +  A,  x  <  —A 

Initialization:  d°  = 

b°  =0,Jfc  =  0, 

ll 

(N 

s 

< 

while  |  Am 

k+l 

u 

=  (KTK  +  Aiyl(KTf  + 

A(dk-bk)), 

dk+x 

=  Sz(uk+x 

+  bk), 

bk+' 

=  bk+uk+ 

l-dk+\ 

Am  = 

-uk+l  -uk , 

> 

k  — > 

k  +  l , 

end 


Algorithm  1  Split  Bregman  algorithm. 


(2) 


3.  Multi-Aerosol  Split  Bregman  Algorithm 

In  this  section  we  generalize  the  basic  split  Bregman  method  described  in  section  2  to 
treat  the  aerosol  unmixing  problem  for  multiple  aerosols  within  a  single  lidar  line-of- 

sight.  Our  starting  point  is  the  MxN  matrix  |G/X.|  that  represents  the  range-resolved 

lidar  data  at  a  single  time-step  after  the  preprocessing  to  remove  the  natural  aerosol 
backscatter  and  transmitter  pulse  effects  as  described  in  [1],  The  row  index  of  G, 
1  <  j  <M ,  labels  the  wavelength,  and  the  column  index,  1  <  k  <  N ,  labels  the  digitized 
range-cell.  We  model  G  as  the  additive  combination  of  the  aerosol  backscatter  from  L 
materials 


Ik  +  }ljk  ' 


1= 1 
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where  pn  are  the  spectral  backscatter  variables  as  a  function  of  wavelength  and  material, 

and  Cik  are  the  concentration  variables  as  a  function  of  material  and  range-cell.  The 
variables  represent  zero-mean  additive  noise.  As  discussed  in  [1],  to  lead  to  an 
identifiable  model  we  must  constrain  p  and  C  for  each  material.  We  adopt  the  choice 
made  in  [1]  of  requiring  that  CIk  >0,  p  jt>  0,  and  the  backscatter  estimates  be  unit 

vectors  for  each  material:  TjPji=L 


Our  goal  is  to  resolve  G  into  L  aerosol  components  using  a  generalization  of  the  Li- 
regularized  split  Bregman  method.  To  that  end  we  perform  the  splitting  by  introducing 
the  constrained  versions  of  p  and  C,  d  and  e,  that  carry  the  Li-nonn  regularizations,  and 
Lagrange  dual  variables  b  and  c  that  enforce  the  constraints  p  —  d  ,  and  C  =  e  within  an 
augmented  Lagrangian  framework.  Setting  the  data  fidelity  term 


H(P'C)  =  \L 


j,k 


“1 2 


/= 1 


we  have  the  augmented  Lagrangian 

L(p,C,d,e,b,c )  =  H(p,C)  +  |c?|  +  ||e||  +  X  (b,p  —  d)  +  Xc  (c,C  —  e) 

2 p  I,  m2  Ar  ||  m2 

H - \p  —  d  H - C-el  , 

2  W  Il2  2  11  1,2 

to  be  minimized  over  p,C,d,e ,  and  maximized  over  b  and  c  such  that  p  >  0,  C  >  0,  and 
||p||,  =  1  for  each  material.  The  parameters  >  0  and  Xc  >  0  control  the  balance 

between  fitting  errors  and  parameter  smoothness,  and  are  set  heuristically  using  a 
combination  of  simulated  and  test  data. 

As  in  the  case  of  the  basic  split  Bregman  algorithm,  the  multi-aerosol  augmented 
Lagrangian  leads  to  a  saddle  point  solution  that  iterates  between  concentration  and 
backscatter  using  three  iterated  steps  for  each  component  given  current  estimates  of  the 
other.  The  steps  for  backscatter  are: 


1 .  pn+x  =  arg  min 


2.  dn+l  =  arg  min 


X, 


ff(p,C")  +  ^(bn,p-d")  +  ^lp-d 


P 


X, , 


\\dl+Xp{bn,pn+l-d)  +  ^-\p^-d 


3.  bn+l  =b"  +pn+l-d'l+ 1 
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Because  of  the  introduction  of  d,  the  problem  splits  into  three  easily  solved  components: 
the  first  step  is  solved  by  calculus,  the  second  step  is  solved  by  shrinkage,  and  the  third 
step  is  explicit. 

The  corresponding  steps  for  concentration  are: 


1. 

2. 

3. 


C  =  arg  min 

c 


n+ 1  • 

e  =  arg  mm 


n+ 1  n  ,  s~i. 

c  —  c  +  C 


n+ 1 


H(pn+l ,0  +  4  (cn ,C-e")  +  ^-\\C-e' 
|4+4(c",0+1-e)  +  ^|0+1-e|[ 


C"+1  ->max(c"+1,0) 


Algorithm  2  gives  the  resulting  unmixing  algorithm.  The  function  SA  is  the  shrinkage 

operator  (2).  The  recursions  in  Algorithm  2  are  initialized  by  prior  estimates  of  the 
backscatter,  if  available,  or  randomly.  Typically,  10  iterations  are  enough  to  give  good 
results. 

4.  Numerical  Examples 

In  this  section  we  illustrate  the  multi-aerosol  unmixing  algorithm  on  simulated  and  actual 
release  data  collected  by  the  US  Army  FAL  (Frequency  Agile  Lidar)  sensor.  The 
calculations  were  done  by  implementing  Algorithm  2  in  MATLAB  on  a  PC.  The 
simulated  FAL  data  were  created  by  injecting  two  overlapping  aerosol  plumes  having 
Gaussian  range-dependence  both  centered  at  1.5  km  with  width  50  m  into  background 
FAL  data  sets.  Those  sets  consisted  of  transmitted  pulse  waveforms  and  received  laser 
backscatter  from  the  natural  atmosphere  as  a  function  of  time-step,  wavelength  index,  and 
range-cell.  The  two  injected  plumes  represent  a  bioaerosol  simulant  (bacillus  BG)  and 
interferent  (kaolin  dust).  The  simulant  plume  was  injected  between  time-steps  200  and 
600,  and  the  interferent  between  time-steps  400  and  800.  A  randomly  varying  peak 
concentration  was  created  from  a  first-order  AR  model.  The  purpose  of  the  simulator 
was  to  generate  data  for  overall  algorithm  perfonnance  verification  including  background 
removal  and  transmitter  pulse  deconvolution.  Because  those  operations  are  not  relevant 
to  the  present  discussion,  we  focus  on  the  comparison  of  the  input  concentration  and 
backscatter  wavefonns  with  the  estimates  from  the  split  Bregman  algorithm  after  data 
preprocessing. 

Figure  2  shows  the  concentration  estimates  from  the  unmixing  algorithm  with 
regularization  parameters  set  at  Ap  =  Ac  =  0.05  .  The  plots  show  the  bioaerosol  estimates 

(left)  and  interferent  (right)  versus  time-step  and  range.  The  calculation  was  initialized 
using  the  mean  of  the  backscatter  spectral  estimates  used  to  train  our  support  vector 
machine  classifier  for  the  bioaerosol  and  interferent  classes. 
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Figure  3  compares  the  input  and  estimated  peak  concentration  over  range  as  a  function  of 
time-step  for  the  bioaerosol  (top)  and  interferent  (bottom).  We  see  that  the  estimates 
track  the  input  throughout  the  run  including  the  plume  overlap  region.  For  comparison, 
Figure  4  plots  the  peak  concentration  estimates  from  a  previous  algorithm  based  on  using 
backscatter  mean  spectra  as  priors  in  a  Bayesian  scheme.  Since  the  prior  as  well  as  data 
densities  were  chosen  to  be  Gaussian,  the  Bayes  approach  in  effect  is  doing  quadratic 
norm  processing.  Although  giving  good  results  when  the  plumes  are  separated  in  time, 
this  processor  fails  to  produce  good  material  separation  in  the  overlap  region. 


Initialize  p°  through  prior  estimates  or  randomly. 

DP  =  1.  Aid  =  0 
while Dp>  s 

e°  =c°  =d°  =b°  =0,n  =  0,  C°  =0,  ||AC||,  =1 
while  ||AC||2  >  £ 

Cn+l  =  ( pnT pn  +  2C/)  *  (gV  +  4  [en  -cn )) , 
Cn+l  -Amax(c”+1,0), 
e"+1  =  SA  (C',+1  +c”) , 
c"+1  =c"  +C"+1  -e"+1  , 

A C  =  Cn+1-Cn  , 
n  — »  n  + 1, 
end 


C  =  Cn+l,  ||A/?||2  =  1,  n  =  0 
while  \\Ap\\  >  £ 


n+ 1 

P 

=  (ccr  + 

V)  (GCT  +  Ap(d" 

n+ 1 

.  1  n+ 1 1  /II 

n+ 1 1| 

P 

I/I 

1^  IL’ 

dn+l 

=  S,(p"' 

+  bn), 

bn+l 

=  b"  +  p"+] 

-dn+x  , 

A p  = 

=  Pn+i~Pn 

n  — > 

n  + 1 , 

end 


P  =  ’DP  =  \P  ~  Pold  ||2  »  Pold  =  P 
end 

Output  p  and  C 

Algorithm  2  Split  Bregman  algorithm  for  multi-aerosol  unmixing. 
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Fig.  2  Concentration  estimates  (bioaerosol  left  and  interferent  right)  from  simulated  data  example. 
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Fig.  3  Split  Bregman  peak  concentration  estimates  for  bioaerosol  (top)  and  interferent  (bottom). 


Figures  5  and  6  compare  the  concentration  range-dependence  (true  and  estimated)  and 
backscatter  wavelength  dependence  at  time  steps  300  and  700  for  the  bioaerosol  and 
interferent,  respectively.  The  results  suggest  that  the  split  Bregman  unmixing  can 
provide  both  low  bias  and  low  variance  estimates. 


Fig.  4  Peak  concentration  estimates  using  a  Bayes  (quadratic  norm)  approach. 
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Fig.  5  Concentration  range-dependence  for  bioaerosol  (top)  and  interferent  (bottom). 
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Fig.  6  Backscatter  wavelength-dependence  for  bioaerosol  (top)  and  interferent  (bottom). 
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As  an  example  of  processing  FAL  data  from  a  mixture  of  aerosols,  we  consider  the 
partially  overlapped  release  TD4065  of  Florida  BG  (a  bioaerosol  simulant)  and  Kuwait 
dust  on  May  31,  2008  in  the  Joint  Ambient  Breeze  Tunnel  (JABT)  at  Dugway  Proving 
Ground,  UT.  The  FAL  sensor  was  located  about  1 .2  km  from  the  entrance  to  the  JABT. 
The  bioaerosol  was  released  between  time-steps  100-476,  and  the  dust  between  times 
261-633.  Data  collected  during  the  first  100  time-steps  were  used  to  estimate  the 
backscatter  from  the  natural  atmosphere.  The  aerosols  were  injected  near  the  far  end  of 
the  tunnel,  and  drawn  toward  the  front  of  the  tunnel  by  large  fans  where  they  were 
exhausted. 

Figure  7  shows  the  concentration  estimates  from  the  unmixing  algorithm  on  these  data 
with  the  bioaerosol  (dust)  plotted  on  the  left  (right).  We  see  the  initial  localized  aerosol 
appearing  at  the  indicated  beginning  of  the  releases.  The  subsequent  time-steps  show  the 
plumes  expanding  to  fill  the  space  between  the  injection  and  extraction  locations.  Figure 
8  plots  the  peak  concentration  over  range  with  the  bioaerosol  (dust)  shown  at  the  top 
(bottom).  We  see  a  clean  separation  of  the  two  aerosol  components  in  these  figures. 

The  estimates  of  aerosol  backscatter  are  plotted  in  Figure  9  as  a  function  of  time-step  and 
wavelength  index  for  the  16  wavelengths  used  in  these  tests.  Up  to  time-step  100  we  see 
only  the  random  unit  vectors  produced  by  the  algorithm  in  the  absence  of  aerosol.  After 
the  bioaerosol  release  (shown  on  the  left)  the  estimates  show  a  consistent  structure, 
although  they  are  rather  noisy  between  time-steps  300-450  because  of  low  concentration 
levels.  The  corresponding  dust  backscatter  estimates  (shown  on  the  right)  also  show  a 
consistent  structure  throughout  the  dust  release. 


Fig.  7  Concentration  estimates  (bioaerosol  left  and  dust  right)  from  release  TD4065. 
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Fig.  8  Peak  concentration  estimates  for  bioaerosol  (top)  and  dust  (bottom)  from  TD4065. 
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Fig.  9  Backscatter  estimates  for  bioaerosol  (left)  and  dust  (right)  for  TD4065. 
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Finally,  the  aerosol  classifier,  a  least-squares  support  vector  machine  (LS-SVM) 
classifier,7  was  applied  to  the  backscatter  estimates.  Our  classifier  consists  of  a  three- 
class  discriminator  for  bioaerosol,  interferent,  and  null  classes.  The  latter  class  was 
included  to  allow  the  classifier  to  perform  a  detection  function  by  rejecting  spectral 
patterns  that  look  like  neither  bioaerosol  nor  interferent.  Three  binary  classifiers  were 
trained  by  the  1-versus-rest  method  for  discriminating  samples  from  each  class  against 
the  pooled  alternatives.  At  each  time-step  the  class  whose  classifier  has  the  highest 
decision  function  value  is  chosen.  Backscatter  estimates  from  both  material  components 
were  processed  this  way,  and  the  results  for  the  decision  functions  are  plotted  in  Figure 
10.  The  results  indicate  that  the  classifier  recognized  both  materials  at  the  correct  times 
during  the  partially  overlapped  releases. 


Fig.  10  LS-SVM  decision  functions  for  bioaerosol  (top)  and  interferent  (bottom)  for  TD4065. 


5.  Summary  and  Conclusions 

The  unmixing  of  elastic  backscatter  data  from  multi-wavelength,  range-resolved  lidar  for 
aerosol  mixtures  has  defied  solution  by  traditional  regularization  methods  using  quadratic 
smoothing  constraints.  A  new  shrinkage-based  L  i  -regularization  method  for  linear 
inverse  problems,  the  split  Bregman  algorithm,  has  been  successfully  applied  to  the 
aerosol  unmixing  problem.  It  was  illustrated  here  on  simulated  and  actual  aerosol 


7  J.  A.  K.  Suykens,  T.Van  Gestel,  J.  De  Brabanter,  B.  De  Moor,  and  J.  Vandwalle,  Least  Squares  Support 
Vector  Machines,  World  Scientific,  2002. 
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mixture  release  data  collected  by  US  Army  personnel  in  field  testing  at  Dugway  proving 
Ground,  UT  using  the  rapidly  tuned  FAL  sensor. 

The  capability  of  resolving  data  from  mixtures  of  aerosols  into  their  backscatter  and 
concentration  components  is  important  to  the  success  of  standoff  sensors  for  biological 
detection  in  realistic  operating  environments  where  interferent  materials  such  as  smoke 
and  dust  will  usually  be  present.  Failure  to  correctly  perform  this  unmixing  can  lead  to 
increased  misclassifications  that  degrade  the  usefulness  of  potential  sensors  employing 
active  detection. 

Other  possible  standoff  sensing  applications  of  this  unmixing  technique  include  (1)  the 
analysis  of  lidar  data  from  chemical  releases  where  the  agent  can  have  both  vapor  and 
aerosol  phases  and  interferent  materials  can  include  byproducts  of  an  explosive  release, 
and  (2)  the  use  of  thermal  imaging  sensors  operating  in  the  long  wave  infrared  spectra 
region. 

Appendix — Augmented  Lagrangian  Derivation  of  Split  Bregman 

Given  the  augmented  Lagrangian  function  L(u,d,b )  in  (1),  we  wish  to  find  a  saddle 
point  ( u*,d*,b  *)  such  that8 

Z(w*,<i*,h)  <  Z(w*,<i*,h*)  <  L(u,d,b*) . 

We  solve  this  saddle  point  problem  by  iterating  between  primal  and  dual  optimizations: 

(uk+1,dk+l)  =  argmmL(u,d,bk)  (primal) 

'  '  u,d  k  / 

bk+1  =  bk  +  uk+l  -  dk+l  (dual) 

For  the  primal  minimizations  holding  the  dual  variable  b  fixed,  we  note  that  because  of 
the  introduction  of  the  d  variables,  the  optimization  over  it  can  be  done  for  fixed  d  by 
differentiation: 


dL(u,d,b)  =K*(Ku_f^  +  Ab  +  A(u_d) 
du 


to  get 


uk+1  =(K*K  +  Aiy1  (K*f  +  A(dk  -bk )) . 


(A-l) 


s  This  sketch  follows  the  split  Bregman  treatment  given  by  Y.  Gui-Bo  and  X.  Xiaohui,  “Split  Bregman 
method  for  large  scale  fused  Lasso,  ”  J.  Comp.  Stat.  &  Data  Analysis,  Vol.  55,  No.  4,  April  2011. 
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Likewise,  although  the  d  variable  optimization  is  still  non-smooth,  it  is  very  easy  to  solve 
by  shrinkage  for  given  u  and  b\ 

dk+l  =  SA(uk+l+bk),  (A-2) 

where  the  shrinkage  function  SA  is  given  by  (2).  Steps  A-l  and  A-2  can  be  iterated  for 
fixed  bk ,  but  we  find  it  more  efficient  (usually)  to  perform  them  only  once. 
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